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We present the interplay between synciironization of unidirectional coupled chaotic nodes with 
heterogeneous delays and the greatest common divisor (GCD) of loops composing the oriented 
graph. In the weak chaos region and for GCD=1 the network is in chaotic zero-lag synchronization, 
whereas for GCD=m > 1 synchronization of m-sublattices emerges. Complete synchronization can 
be achieved when all chaotic nodes are influenced by an identical set of delays and in particular for 
the limiting case of homogeneous delays. Results are supported by simulations of chaotic systems, 
self-consistent and mixing arguments, as well as analytical solutions of Bernoulli maps. 

I. INTRODUCTION 

Synchronization, complex networks, and chaotic dynamics with delay couplings are emerging phenomena, and 
concepts which have fascinated scientists for decades. These phenomena are ubiquitous in nature and play a key 
role in almost all fields of science including biology, ecology, physics, climatology, sociology and technology [ll-yl- 
However, each one originates and is governed by different features and rules. For instance, the description and 
classification of complex networks are often based on their statistical properties, such as degree distribution, average 
degree and degree correlations The observation that real networks have degree distributions that are very 

different from those of classical random graphs was the starting point for the recent explosion of interest in complex 
networks 0, Q . By contrast the dynamics of processes defined on networks are closely related to the spectrum of 
an appropriate connection operator; a prototypical example is chaos synchronization Q, which crucially depends on 
the extreme eigenvalues of the graph Laplacian [sl-fioj. Is there an interplay between the statistical properties of a 
network and its extreme spectral properties? Over the last few years a number of papers have reported correlations 
between the synchronization of a network and its degree of homogeneity [TTI - [l3j , clustering coefficients [ill , degree 
correlations fig , average degree and degree distribution [l6j . The literature even reports some conflicting trends, 
e.g. synchronization is amplified/damaged by increasing the degree of homogeneity, or adding a few shortcut links 
enhances/reduces the level of synchrony [ill Ull, [3 . The effects of distributed time delays and coupling strengths in 
neurons and coupled oscillators have been studied, showing that signal transmission can be seriously influenced by 
the distributed time delays and significant delays may result in synchronization [l7l - [20| . However, a general theory 
which relates the structure of the network to its synchronization properties, especially in the case of networks with 
heterogeneous time-delay couplings and where the dynamics of the network is chaotic remains elusive. 

In this Letter we report that synchronization of chaotic networks composed of identical nonlinear units with het- 
erogeneous unidirectional time-delayed couplings is governed by the GCD of the length of directed loops of the 
corresponding graphs. We consider strongly connected oriented graphs, i.e. cases where there is a path from each 
node to every other node in the graph. Each directed edge of the graph corresponds to a coupling with delay time 
fc^r, where ki is an integer and r is the time unit of the coupling delay. The length of a directed loop is the sum of 
all integers ki along the loop, and the GCD of lengths of all loops determines the dynamic properties of the network. 
Two main types of chaotic synchronization can be achieved as a function of the value of the GCD in the limit of weak 
chaos, i.e. a small positive maximal Lyapunov exponent Xmax- For GCD=1, zero-lag synchronization (ZLS) among 
all nodes of the oriented graph can be achieved, whereas for GCD=m > 1, nodes are partitioned into m-sublattices 
(clusters), where all nodes belonging to a sublattice are in ZLS. This synchronization state is different from Chimera 
states, where a network splits only into synchronized and desynchronized sub-populations [2ll.[2^. 



II. GCD-SUBLATTICES 



The results of heterogeneous oriented small networks exemplifying the role of the GCD are presented in Fig. 1 for 
the chaotic Bernoulli maps (BM) [2^ , as well as for the Lang-Kobayashi equations (LKE) which are a good model for 
the intensity dynamics of coupled chaotic semiconductor lasers [2^ and are explicitly given in references [2^, [2^ . 
The strength of a directed coupling is denoted by a and e for the LKE and BM, respectively, and for simplicity we 
assume that the sum of incoming equal strength couplings to each unit is identical, i.e. for the BM case and for a 
unit with q incoming couplings the dynamics is xl^i = (1 — ()f{x\) + ^/q'Yl'j fi^l-k r)^ where x] is the state of node 
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FIG. 1: Synchronization of homogeneous/heterogeneous oriented chaotic motifs. Nodes with tlie same color are synchronized 
and the number of colors represents the number of SLs. Small colored nodes in heterogeneous motifs are presented for illustration 
only, (a) Homogeneous motif consists of two connected loops of sizes 9r and 12r with 3-SLs. (b) Heterogeneous motif similar 
to (a) with 2-SLs. (c) Similar to (b) with 3-SLs. (d) Homogeneous motif similar to (a) in ZLS. (e) Heterogeneous motif in 
ZLS, where tick marks indicate unit distances r. (f) Three connected heterogeneous loops of sizes 6r, 15r and lOr, where 
GCD(6,10,15) =1 and the network is in ZLS. (g) Pairs of connected loops taken from the three loops of (f), where each pair is 
in SLs in contrast to ZLS in (g). 



i at time t and f{x) = {ax)modl which is chaotic for a > 1. Throughout this report r is equal to 40 steps and 10 ns 
for BM and LKE , respectively, unless otherwise indicated. 

Figure la depicts a network composed of two connected loops of lengths 9r, 12r and sixteen nodes. Since 
GCD(9,12)=3, a chaotic synchronization consisting of 3 sublattices (3-SLs) is expected, and is depicted in Fig. la 
where each SL is represented by a different color. Solving the master stability function Q for the sixteen BM reveals 
that for a=1.02, for instance, and a coupling strength e > 0.47 the 3-SLs is a stable solution. Figure lb presents two 
loops with identical total delays as in Fig. la, but with only 5 nodes, where the heterogeneous square/triangle consists 
of (4t, r, 4r, 3t) / (4r, 2r, 3r) delays between connected nodes. Results of simulations for BM with a=1.02 and e > 0.79 
as well as for the LKE equations with p = 1.01 and a = 17ns~^ indicate a 3-SLs solution, where the "color" of the 
5 nodes are fixed by the corresponding nodes in Fig. la. Figure Ic presents a similar heterogeneous square/triangle 
consisting of (4r, 2r, 3r, 3r)/(4T, 2t, 3t) delays between directed connecting nodes. Now the 5 nodes form only 2-SLs 
in line with the corresponding nodes of Fig. la. 

The main difference between homogeneous and heterogeneous oriented networks is the level of synchronization 
between nodes belonging to the same SL. For homogeneous networks ki ~ 1, Figs, la and Id, complete synchro- 
nization is achieved such that nodes belonging to the same SL have an identical chaotic trajectory. By contrast, 
for heterogeneous networks, complete synchronization is not a solution of the dynamics, since units belonging to the 
same SL can be connected to preceding units by different delays and chaotic signals are not periodic. Nevertheless, in 
the weak chaos region, a solution with a remarkable ZLS between trajectories of nodes belonging to the same SL of 
the heterogeneous network emerges, where the level of synchronization is enhanced toward a complete ZLS as chaos 
becomes weaker. For the parameters of the heterogeneous networks synchronization measured by cross correlation 
[26| at zero time shift C ^ 0.9, where for BM it was averaged over a window of 10"^ steps and for LKE it was averaged 
over windows of 10 ns, excluding the low frequency fluctuations regions [1^ [2^. Note that correlation below 0.5 is 
the typical outcome of in vivo experiments in neuroscience, for instance, when measuring synchronization of activity 
among different areas of the brain [27 1. 

The ZLS for a homogeneous network consisting of the two connecting loops of delays 3t and 4r (the same ratio 
between the two loops, 9r and 12r, of Figs, la-c) is presented in Fig. Id. By rescaling delays with the minimal 
delay in the network, 3t, the effective relative sizes of the loops are 3 and 4 and GCD(3,4)=1. In fact, simulations of 
the network with BM as well as the analytical solution of the master stability function indicate a complete ZLS for 
a=1.02, for instance, and e > 0.12, see Appendix. Similar results were obtained in simulations for LKE with p = 1.02 
and a = 13ns~^. ZLS is also achievable for heterogeneous networks as exemplified in Fig. le. The network consists 
of 5 nodes and loops of sizes 5,7,8,9, 10 and 12 with GCD=1. Simulation results indicate that C is at least 0.9 for 
BM with a=1.02 and e > 0.89 and for LKE with p = 1.02 and a = 17.5ns-^. 

A more complex heterogeneous network consisting of three connecting loops for a total delay of 6r(T, 2r, 3t), 
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10T(r,2r,3T,4T) and 15T(r, 2r, 2r, 4r, 3r, 3r) is depicted in Fig. If. The GCD(6, 10,15)=! and the entire network 
was expected to be in ZLS, as was confirmed in simulations for BM with a=1.02 and e = 0.92, for instance, where 
C ~ 0.9. Figure Ig depicts the synchronization of networks consisting of only a pair of connecting loops of Fig. If, 
(6,15), (6,10) and (10,15), where the chaotic behavior is GCD(6,15)=3-SLs , GCD(6,10)=2-SLs and GCD(10,15)=5- 
SLs, respectively, as was confirmed in simulations for BM. These results indicate that each motif does not maintain 
its behavior in the entire network, which thus cannot be simply described as a "Lego" of connecting components 
with given chaotic modes of activity. Note that motifs can be connected either by common delay couplings or by a 
single node. Hence the role of GCD is a global decision which in general cannot be deduced from local topological 
or geometric properties of the network. An exception is the case where GCD=1 for two local loops, which is enough 
information to deduce that ZLS takes over the entire chaotic behavior of the network. Nevertheless, the GCD in 
general can induce long-range effects such that addition, deletion, changes in delays of existing couplings can affect 
the entire chaotic state of the network and in particular correlations of remote nodes. 



III. CHAIN AMPLIFICATION 



Synchronization of homogeneous networks is achievable for weaker chaos, e.g. a smaller a for BM, as the size of 
loops increases. This trend can be attributed to the typical emergence of longer chains in a network composed of larger 
connected loops, where the largest Lyapunov exponent (amplification) scales in the first order approximation linearly 
with the number of nodes that constitute the chain, see an illustration in Fig. 2. Hence the emergence of longer 
chains requires a weaker chaos to maintain complete synchronization. A similar trend is applicable for heterogeneous 
networks, which resemble homogeneous networks with additional intermediate units which result in longer chains. In 
addition, to maintain a high level of synchronization for the heterogeneous case a weaker chaos is required than for 
the corresponding homogeneous network characterized by complete synchronization. Both these trends were observed 
in simulations and for some heterogeneous networks the master stability functions were explicitly solved. Hence 
synchronization of chaotic networks is expected in general to be found in the weak chaos region. 




FIG. 2: Two Bernoulli chains with the same length, 4r: The diflerence between the state of the last unit is compared between 
two close initial conditions, x and a; + A, for the first unit. For e — >■ 1, for instance, the internal dynamics is negligible and the 
difference for (a) is a* A whereas for (b) it is aA. Hence, the Lyapunov exponent of the entire chain is ~ 41og(a) and ~ log(a) 
for (a) and (b), respectively. For directed graphs composed of large loops, longer chains are more likely to appear and in order 
to compensate the amplification of such chains, compare to graph with smaller loops, a smaller slope a has to be selected. 



IV. MULTIPLE DELAYS 



For networks with GCD=1, complete ZLS is achievable either for the case where all nodes are influenced by an 
identical set of heterogeneous delays or for homogeneous delays, ki = 1. Figure 3a depicts the heterogeneous motif 
with minimal number of nodes with complete ZLS. This motif consists of 3 nodes, each one of which is connected to 
the preceding node with two delays, t and 2r and the motif consists of loops of 3r, 4r, 5t and 6t such that GCD=1. 
Similarly, ZLS can be found where 2t is replaced by 3r, 5r etc., where GCD=:1. Note that the topology of this motif 
is similar to Fig. le; however, the geometry is different. Figure 3b depicts the homogeneous motif with minimal 
number of nodes, which consists of 4 nodes composed of loops of 3r and 4t such that GCD=1. Results of simulations 
for BM as well as the analytical solution of the master stability function for both networks indicate ZLS for a=1.02, 
e > 0.1 for Fig. 3a and e > 0.12 for Fig. 3b. Similarly, correlation close to oiie (~ 0.99) in the ZLS state was measured 
in simulations out of the low frequency fluctuations regions for LKE [2^, [2^ for p=1.02 and a = 147is^^ for Fig. 
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3a and a = 17.5ns~^ for Fig. 3b. Note that complete ZLS can be found in a smaller heterogeneous motif than in 
homogeneous motifs and furthermore multiple delays enlarge the region, e.g. e in Bernoulli, where complete ZLS is 
achieved. 




(b) 




FIG. 3: Oriented motifs with minimal number of units exhibit complete ZLS. (a) Heterogeneous motif consists of 3 nodes, 
where tick marks indicate unit distances, (b) Homogeneous motif consists of 4 nodes. 



MIXING ARGUMENT 



The role of the GCD can be best understood by the self-consistent argument that a necessary condition for a chaotic 
synchronization is that each node is driven by the same set of " colors" , where for heterogeneous networks the missing 
nodes arc artificially inserted, similar to Figs, lb, Ic and Ig. The trivial solution is always one "color", ZLS; however, 
the alternative solution consists of exactly GCD "colors", GCD-SLs. An attempt to consistently color the network 
with a smaller number fails, since a contradiction emerges where nodes with the same color have different drives. As 
for bidirectional chaotic networks |23|, the results always indicate that with a lack of self-couplings, SL always takes 
over ZLS when it is a consistent solution. Another interesting argument that accounts for the emergence of ZLS in 
oriented graphs is the mixing argument which was first proposed for bidirectional networks with multiple delays (28j 
and was recently observed in an experiment on mutually coupled chaotic lasers (29j . Figure 4a depicts the adjacency 
matrix, G, of Fig. 3b. The iOth power of the matrix G indicates that it is a primitive matrix where each one of 
the four nodes receives at time t an input from all four nodes, including the node itself, from time step t — 40, where 
time steps are normalized with r. Furthermore the drives for all nodes are identical indicating that only ZLS is a 
consistent solution. Figure 4b depicts a homogeneous graph consists of two connected loops of lengths 3r and 6t with 
3-SLs and its adjacency matrix. The AOth power of the matrix G indicates 3-SLs: rows/colunms (1,4,7) are identical 
as well as rows/columns (3,5) and (2,6). All nodes belonging to the same clusters are mixing the same information 
from t — 40. The structure of the matrix, G^°, also indicates that for instance the cluster (1,4,7) is driven by the 
cluster (3,5). 

A change in only one directed coupling to bidirectional (mutual) has a dramatic effect on the synchronization 
pattern of a heterogeneous/homogeneous network, since a loop of size 2 is now embedded in the network. As a 
result, in the case where 2 is a common divisor of all loops of the network 2-SLs takes over, otherwise ZLS is the 
solution. The effect of one bidirectional coupling is exemplified in Fig. 5 for homogeneous networks, where similar 
results were obtained for heterogeneous networks as well. Figure 5a consists of two connected directed loops of 8t and 
12t with GCD=4. After one directed coupling is converted to bidirectional (red coupling) 2-SLs is the only possible 
synchronization, e.g. a=:1.01 and e > 0.24. Figure 5b consists of two connected loops 6t and 9r with GCD=3, where 
a bidirectional coupling drives the network from 3-SLs to ZLS, e.g. a=1.01 and e > 0.04. Note that the effect of a 
bidirectional delay, r, is similar to the same network with additional self-coupling 2r, as well as the reverse operation 
where the loop of length 6t in Fig. 5a, for instance, can be replaced by self-coupling of 6t, which were confirmed in 
simulations. 
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FIG. 4: (a) The adjacency matrix, G, for the network of Fig. 3b and its 40th power (first two leading digits of each matrix 
element), (b) Two connected loops of 3t and 6t with 3-SLs. The adjacency matrix, G, and its iOth power. 



VI. ANALYTICAL RESULTS t ^ oo 



For homogeneous oriented Bernoulli networks in the limit of infinite delays, r — > ao, the role of the GCD can be 
established analytically. The corresponding equations for the BM are 

xl = (l ~ e)fixU) + ^Y. Gvfi4-r) ■ 

The adjacency matrix dj represents the couplings and their weights in the oriented network and we assume that the 
sum of the incoming signals to each unit is equal to one, Gij = 1 . These special types of non- negative matrices are 
known as stochastic matrices and play a central role in Markov chain processes, where many of their mathematical 
properties are known [30| . 

For Bernoulli networks with homogeneous infinite delay couplings the master stability function depends solely on 
the eigenvalue spectrum of the adjacency matrix G [1^. More precisely, the matrix G always has an eigenvalue 70 = 1, 
which determines the Lyapunov exponent tangential to the synchronization manifold (SM) and does not affect the 
stability of the synchronization. The stability of the SM is given by 



7 < e-^— " 
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FIG. 5: Change of one directed coupling to bidirectional, (a) Two connected loops of 8r and 12r, where one bidirectional 
coupling (red) changes synchronization from 4-SLs to 2-SLs. (b) Two connected loops of 6r and 9r , where one bidirectional 
coupling (red) changes the synchronization from 3-SLs to ZLS. 



where 7 is the second largest modulus of the eigenvalues of G, and Xmax is the largest Lyapunov exponent. Hence, 
a sufficient condition for the stability of the SM is determined by a non-zero eigenvalue gap (7 < 1). Markov chain 
theory now makes it possible to derive mathematical statements such as: (a) bi-partitc networks as well as directed 
rings have 7 = 1 and complete ZLS is unstable; (b) a network where each node is connected to any other node and 
GCD=1, 7 < 1, hence complete ZLS is possible; (c) for a similar network with GCD=m, G™ has a block structure of 
m blocks and the network is in m-SLs. 



VII. CONCLUDING REMARKS 



The interplay between GCD and types of synchronization was found to be robust for chaotic networks with some 
heterogeneity of coupling strengths as well as for networks with self-couplings which is nothing else but additional 
directed loops. Since the GCD is typically a global quantity, synchronization of networks cannot be simply described as 
a "Lego" of connecting components with given synchronization modes; hence, this casts some doubt on the importance 
of their statistical properties [sij. Small changes in geometry and topology can dramatically alter the number of GCD- 
SLs, such that addition/deletion of a precise coupling can serve as a remote switching mechanism in the network. 
However, it should be borne in mind that the sparseness of networks is crucial for the richness of synchronization 
modes, since for highly dense networks the GCD is expected to be typically one. 



VIII. APPENDIX 



We exemplified calculations of the master stability for the network depicted in Fig. Id. The dynamical equations 
are given by 

< = (l-e)/«_i)+e^G,,/(.T^„_J, » = !,..., 5 (1) 
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(2) 



For a small perturbation 5x1^ from the trajectory x'j , one can linearize the equations 

5< = (l-e)a<5<_i+e^G„afe^„_„ z=l,...,5 (3) 

j 

and using the ansatz 6x1^ = c^Sxq one finds 

cSxq = {1 — e)a5xQ + e'^^Gijac^~'^Sx-'Q, i — l,...,5. (4) 
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The matrix G has five eigenvalues, 7^. Substituting the eigenvalues into this equation we get the following characteristic 
polynomial: 

- (1 - e)ac^-i - ea7,; = 0, i = l,...,5 (5) 

where c = expA + i(/) and A = ln\c\ is the Lyapunov exponent. The polynomial determines the entire spectrum of 
the Lyapunov exponents of the system. However, only the exponents transversal to the synchronization manifold are 
important for the stability of the synchronization. The polynomial is of order r, therefore it has r solutions, Lyapunov 
exponents, for each eigenvalue, 7^. The synchronization is stable when 7^ < Vi except for those correspondents to 
71 = 1 , which is parallel to the synchronization manifold. Hence checking only the maximal one for each eigenvalue 
is sufficient. We solved the polynomial using Matlab for a = 1.02 and r = 40 and found that the synchronization is 
stable, i.e. for each of the four eigenvalues the maximal Lyapunov exponent is negative, for e > 0.12. 

A discussion with Fabian Wirth is acknowledged. The work of LK. is partially support by the Israel Science 
Foundation. 
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